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In systems removed from equilibrium, intrinsic microscopic fluctuations become correlated over 
distances comparable to the characteristic macroscopic length over which the external constraint is 
. exerted. In order to investigate this phenomenon, we construct a microscopic model with simple 

^ ' stochastic dynamics using lattice gas automaton rules that satisfy local detailed balance. Because of 

the simplicity of the automaton dynamics, analytical theory can be developed to describe the space 
and time evolution of the density fluctuations. The exact equations for the pair correlations are 
solved explicitly in the hydrodynamic limit. In this limit, we rigorously derive the results obtained 
phenomenologically by fluctuating hydrodynamics. In particular, the spatial algebraic decay of the 
^ ■ equal-time fluctuation correlations predicted by this theory is found to be in excellent agreement with 

'^J ' the results of our lattice gas automaton simulations for two different types of boundary conditions. 

2^ . Long-range correlations of the type described here appear generically in dynamical systems that 

' exhibit large scale anisotropy and lack detailed balance. 

^ . PACS numbers: 05.20.-y, 05.40. -fj, 51.10.+y 

^ ■ I. INTRODUCTION 

O I In a hydrodynamic system under non-equilibrium conditions the fluctuations of the densities of conserved quantities 
are correlated over large distances, as confirmed by recent experiments performed by Law, Sengers et al The 
Q ' presence of long-range correlations in systems removed from equilibrium had been predicted by kinetic theory , by 
^ , non-equilibrium statistical mechanics ||^, and by phenomenological theories, such as fluctuating hydrodynamics [Q. 
These correlations decay algebraically over distances comparable to the size of the system. They appear generically in 
systems subject to non-equilibrium constraints PJ^Il^ and are a consequence of the existence of conserved quantities, 
the absence of detailed balance, and the presence of spatial anisotropy. The long-range nature of these correlations is 
remarkable in as much as it is manifested in fluids where the interactions are short-ranged, and under conditions for 
which the fluid is far from critical points or hydrodynamic instabilities. The fact that equilibrium correlations remain 
short-ranged away from critical points is a consequence of precise cancelations of the effects of noise sources Q . As 
soon as this balance, which is characteristic of the equilibrium state, is lost the correlations may become long-ranged, 
as they generically do. 

Our objective in the present work is to put forth and analyze a simple microscopic model, which nonetheless possesses 
sufficiently complex dynamics to exhibit this type of long-range correlations. In particular we construct a lattice gas 
automaton (LGA) corresponding to a collection of random walkers. The "particles" move on a regular array at discrete 
time intervals and interact via an exclusion principle, a constraint which acts as a sort of hard-core potential in the 
lattice. Particles also enter collisions whose outcome, while conserving the number of particles, is otherwise entirely 
random. These collision rules satisfy a local detailed balance relation, and conserve momentum globally, in a statistical 
way, but not locally. At a global scale, detailed balance is absent because of the imposition of non-equilibrium 
constraints. The automaton dynamics naturally lends itself to a hierarchical description [^7|-^. Microscopically, 
particles propagate between adjacent nodes and experience collisions. From a "macroscopic" point of view, the 
evolution of the automaton is given by a diffusion equation. An intermediate "mesoscopic" description accounts for 
the statistical properties of the fluctuations, which correspond to deviations from the average behavior arising from 
the microdynamics. In equilibrium, the fluctuation correlations are localized on a single node. In the presence of a 
density gradient, which maintains the automaton away from equilibrium, we can distinguish two contributions in the 
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correlation function of the particle number fluctuations: A local equilibrium one, which is short-ranged (in our model, 
it is strictly localized on one lattice node), and a long-range term, which decays algebraically with a characteristic 
length on the order of the system size. The simplicity of the dynamics of this automaton at the microscopic level 
makes it possible to develop an analytic description not only for the evolution of the densities of conserved quantities 
of the automaton (in this case, solely the number of particles per node), but also for their fluctuations. Thus, we 
derive exact equations for the evolution of pair correlations and solve them in the hydrodynamic limit. The lowest 
order approximation in a perturbative scheme using the inverse of the size of the system as an expansion parameter, 
leads to the same expressions for pair correlations as the phenomenological theory of fluctuating hydrodynamics. The 
validity of this approximation is tested against the results of simulations of the dynamics of the automaton. It is seen 
that even for very small automata (that is, even rather far from the thermodynamic limit) fluctuating hydrodynamics 
provides a very accurate description of the statistical properties of the fluctuations. 

The present model system is similar to a cellular automaton proposed by Kawasaki |2^] and studied by Spohn , 
who also derived the existence long-range correlations from the automaton microscopic dynamics. From a mesoscopic 
viewpoint, both the cellular automaton studied by Spohn and the lattice gas automaton (LGA) investigated here 
represent two different microscopic realizations of a stochastic equation analogous to that studied by Nicolis and 
Malek-Mansour and Garcia et al ||2^ to describe heat transport in a rigid conductor subject to a temperature 
gradient. In this case, the conserved quantity is the energy density (as measured by the local temperature), and the 
temperature fluctuations are described by a Fourier (heat diffusion) equation, to which a random heat current with 
a local equilibrium form is added. 

Long-range correlations of the kind described in this paper are also present in lattice gas automata (LGA) with 
collision rules that violate detailed balance. These automata attain a homogeneous equilibrium state which is non- 
Gibbsian. Initially automata violating detailed balance were constructed to simulate hydrodynamics at high Reynolds 
numbers . It was realized later on that these models are intrinsically interesting as a paradigm for driven systems 
1^,^, and that they exhibit all the wealth of behavior characteristic of systems removed from equilibrium. In 
particular, they also exhibit algebraically decaying correlations, which have been studied in great detail [ p"2|jTg|j25| ] . 

Our work complements that carried out by these authors and provides a systematic comparison between the 
theoretical description and simulations. One of the main objectives of this paper is to derive the exact equations for 
the evolution of the hydrodynamic variables and their fluctuations and to establish the connection with fluctuating 
hydrodynamics. Once this program is realized, we can describe in detail the mechanism by which the long-range 
correlations are built up from the local microdynamics. We also discuss the validity of the local-equilibrium hypothesis, 
a basic assumption (usually justified a posteriori) in the theory of fluctuating hydrodynamics. 

In section II, we review the problem of heat transport in a rigid conductor using the phenomenological approach 
of fluctuating hydrodynamics. The novelty of the present treatment is that the postulated Landau equation includes 
explicitly the heat reservoirs that maintain the temperature gradient across the conductor. In this scheme it is 
possible to discuss the effect of boundary conditions in the stochastic equation rigorously. Given the long-range 
nature of the correlations, boundary effects should be non-trivial. It is argued that the implemented boundary 
conditions for the automaton described in section III (vanishing of long-range correlations at the boundaries of the 
system) correspond to the paradigm of a system in a quasi-stationary steady state maintained by contact (diffusive or 
thermal) with reservoirs. Section III constitutes the main body of the paper. We construct a two-dimensional lattice 
gas automaton whose collision rules satisfy a local detailed balance relation. Global detailed balance is broken by 
imposing non-equilibrium constraints. The equations for the evolution of the average number of particles per node 
and for the corresponding fluctuations are derived from the microscopic propagation and collision rules. Simulations 
of the dynamics of the automaton demonstrate the accuracy of the theoretical description. Section IV contains a 
summary of results and some concluding comments. 



II. FLUCTUATING HYDRODYNAMICS. 



We consider the problem of heat transport in a one-dimensional rigid conductor. Our starting point is the isolated 
system analyzed by Procaccia et al in Appendix B of Ref. (see Fig. 1). We proceed in successive steps: First, 
the heat diffusion equation for the system depicted in Fig. 1 is solved. Then, we discuss the conditions under which 
such system supports a quasi-stationary non-equilibrium steady state (quasi-NESS) with a linear temperature profile 
for the central portion. Finally, we solve the stochastic equation, which is constructed by adding a random heat flux 
with local equilibrium form to the Fourier equation, in the limit in which the quasi-NESS is established. 
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FIG. 1. Quasi-NESS temperature profile. 
Consider the Fourier equation in one dimension, 



T{x,t) 



(1) 



where k{x) — Ki {9{—x) + 9{x — L)) + K2 0{x) 6{L — x) is the piecewise constant thermal diffusivity ( ^(a;) is the 
Heaviside step function). The initial condition is 



T(x,0) = Tie{-x) + T2 6{x)e{L-x) +Tse{x-L), 



(2) 



where T2 = = Ti + /3-| , and (3 = ^^^^ is the magnitude of the temperature gradient, and a measure of 

how far from equilibrium the system is. 

Rather than directly solving Eq. (1), we consider the equivalent equation for the spatial derivative of the temperature 
profile 



dt dx 

with the corresponding initial condition 

dT{x, 0) 
dx 

The formal solution of Eq. (3) is given by 

dT{x,t) 
dx 
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dT(x, 0) 
dx 



The following step is to find the eigenfunctions Pk{x) of the operator ^^k{x) : 

^ k{x) Pk{x) = - e Pk{x), 

which are given explicitly in Appendix A. Here, we just quote the result for Pk{x) in the limit e = 
for the case of semi-infinite reservoirs: 

Pk{x) « — cos^e{x)6{L-x) + 0{e), 



and 



4 ~ ^ ^{k 



(3) 
(4) 
(5) 

(6) 
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when e — > 0. 

Using the spectral decomposition of the identity 

S{x-xq) = 



dk kIxq) Pfe(xo) Pkix), 



we can give an exphcit solution for Eq. (3) with the initial conditions given by Eq. (4) 

dT{x,t) ^ L f ^2 , 1 , ^, 

—^^P - e^pi^—K{x)tj [S{x-L)+d{x)] 



(9) 
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This solution is valid for t ^ -k — and e 
We note that for an intermediate time regime 
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the derivative of the temperature profile is well approximated by 

^ Tix.t) « p9{x) e{L~x), 



(11) 



which corresponds to the long-lived quasi- NESS depicted in Fig. 1, previously derived in Rcf. p^ . 

In order to understand what are the proper boundary conditions in a non-equilibrium system when fluctuations are 
considered we now proceed to solve the stochastic equation 



d_ 

dx 



g{x,t), 



(12) 



obtained by adding to the deterministic heat-diffusion equation (Eq. (1)) the gradient of a random heat flux g{x, t), 
which is assumed to be Gaussian white noise with the properties: 



{gix,t)) =0, 

gix,t)gix,t)) ^2^k{x) {T{x,t)f 5{x - x ) 6{t - t ), 



(13) 



where ks is the Boltzmann constant, Cp the heat capacity per unit volume, and {T{x,t) ) is the temperature 
profile which is a solution of the deterministic equation, Eq. (1). This Landau equation (Eq. (12)) is constructed 
phenomenologically by requiring that at homogeneous equilibrium {{T{x,t) ) = T^q) the pair correlations are given 
correctly. The extension to a non-equilibrium situation makes use of the hypothesis of local equilibrium. It is assumed 
that the thermodynamic variables that characterize the system are well defined locally. The non-equilibrium steady 
state is thus viewed as a state where these thermodynamic variables vary slowly in space (on a hydrodynamic scale, 
which is much larger than the scale at which the underlying microdynamics take place). In general, the derivation of 
Eq. (12), even in an approximate way, is rather complicated. In the following section we analyze a diffusive lattice 
gas automaton for which we derive from the actual microscopic dynamics a stochastic Landau equation. 
The equal time pair correlation function of the temperature fluctuations is defined as 

C{x,x';t) ^ (5T{x,t) 5T{x ,t)^ , (14) 

with 6T{x,t) = T{x,t) — {T{x,t)). This correlation function is the solution of the differential equation 

_C{x^,,t)^[^-n{x)- + g^<x)g^) C{x,x;t) 

+ ^ ^ ^ <x) {TMf Six-x). (15) 
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This equation is obtained by solving formally Eq. (12), constructing the equal time correlation fimction, and using 
(13) to evaluate one of the time integrals; differentiation of the resulting expression with respect to time yields Eq. 
(15). 

There are two contributions to the correlations: a local-equilibrium contribution (noted LE), and the remainder 
(noted LR), which is long-ranged 

C{x,x;t) = C^^{x,x',t) + C^^{x,x',t), (16) 
The local equilibrium contribution has the form 

C^''{x,x;t) = ^{T{x,t)f 6{x-x'), (17) 

and the long-ranged term obeys the differential equation 

Equation (18) corresponds to a diffusion equation in two dimensions with a source term on the line x = x , which is 
proportional to the square of spatial derivative of the temperature profile. Assuming that the system is initially in a 
local equilibrium state ( C^^{x,x ) = ), the formal solution of Eq. (18) is 



X I J S{x-x). 



(19) 



The explicit solution is obtained by finding the spectral decomposition of the operator §^k{x)^ 

Qk{x) = - Qk{x). (20) 

ox ox 

The details of the derivation are given in Appendix A. Here, we just quote the result for semi-infinite reservoirs in 
the limit e = .f^ ^ 

Qk{x) « sm^0{x)e{L-x), (21) 

with 

oo , 

bl^^ E Sik-n^-^) (22) 

n=—oo 

Hence, in the limit e^O the operator ^'^(2;)^^ has the same spectral decomposition as the operator ^2^^ with 
Dirichlet boundary conditions at a; = 0, L. Assuming the system is in the quasi-NESS described at the beginning 
of this section (i.e. {T{x,t)) = Ts{x) = Ti e{-x) + {Ti + Px) e{x) 6(1 - x) + T3 61(0; - L) ), the long-range 
contribution to the pair correlation function is 

C-(,. ;„ = ji, ^1 f -1 (1 - ,23, 

^ n=l ^ ' 

in the limit e = \/~ ^ and for times t <C t-j— i — • In the time regime during which the quasi-NESS is 



established ( 2^2^^ ^ t -C ^27?^) the long-range steady state correlations are of the form 



5 



^ n— 1 
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kB 0_ 

Cp L 



^^{L~x') X e{x-x') + [L-x) X 0{x - , (24) 



to lowest order in e. The correlations are proportional to the square of the temperature gradient and they decay 
linearly over a distance comparable to the system size. 

It is important to stress at this point that different boundary conditions from the ones described would give rise to 
long-range correlations of different form. However, if we impose the gradient by contact with macroscopic reservoirs, 
the appropriate boundary conditions arc that the fluctuations have a local equilibrium form at the edges of the system. 



III. DIFFUSIVE LATTICE GAS AUTOMATON. 



We construct a model system with stochastic microdynamics which, at the mesoscopic level, can be approximately 
described by the phenomenological theory presented in the previous section. The model belongs to the class of lattice 
gas automata, where particles move on a regular array at discrete time intervals. Since we are interested in the 
transport of a scalar quantity, it is sufficient in two dimensions to consider the case of a square lattice to obtain a 
diffusion equation at the macroscopic level. 

The system analyzed is composed of particles moving on a two dimensional square lattice C of dimensions 
(Lj. + 1) X Lj, |2^. There are four channels per lattice node. Each channel corresponds to the particle velocity 
pointing towards one of the four nearest neighbors ( i = {right) ^ 1 {up), 2 {left), 3 {down) ). Channel i at 
node r G C is occupied if there is a particle at r with velocity Cj. At time t, the state of the automaton is thus 
described by a set of bits 

{n,(r,t); r G £, i = 0,1,2,3}, (25) 

where ni{r,t) is equal to 1 (0) if channel i of node r is occupied (empty). There is an exclusion principle, which 
is the equivalent in the lattice gas of a hard-core potential preventing two particles from simultaneously occupying 
the same channel of a given node. 

The microscopic evolution of the automaton takes place in two stages. There is a propagation step, in which each 
particle in the automaton moves to the same channel of a neighboring node according to its velocity (see Fig. 2(a) 
). For instance, if there is a particle in the channel "up" of node r, the propagation step will move it to the channel 
"up" of the node r -t- y, where y is a unit vector in the Y direction. 

The second stage of the dynamics can be viewed as a "collision" step consisting of a random redistribution of the 
particles at each node (see Fig. 2(b) ) such that configurations compatible with the values of the conserved quantities 
(in our case, the particle number only) are equally probable outcomes of the collision step. This choice for the collision 
rules greatly simplifies the subsequent analytical derivation of the equations describing the macro- and mesoscopic 
evolution of the automaton. 



FIG. 2. (a) Propagation step in a portion of the automaton. The motion takes place on a square lattice with spacings 
Aa; = Ay = 1. A dot on a link, close to node r, indicates that there is a particle occupying that node with a unit velocity along 
the direction indicated by the link. 
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FIG. 2. (b) Collision step. The in state with m = 0, ...,4 on a node yields any of the corresponding states with the 
same number of particles with probability ^7^, where u^g-y is the multiplicity of the state (i.e. the number of compatible 

configurations). 



Prom this set of rules that govern the dynamics of the automaton we can derive a hierarchy of contracted descriptions 
containing progressively smaller amounts of information about the details of the "microscopic" configuration of the 
automaton. The coarsest description is obtained by averaging the equation for the evolution of the channel occupation 
number over the stochastic dynamics. For the class of automata investigated this procedure yields an exact Boltzmann 
equation, which, in the hydrodynamic limit leads to a diffusion equation. Section III. A contains a detailed discussion 
of this "macroscopic" level of description. A more refined (mesoscopic) description is obtained by incorporating the 
statistical properties of fluctuations. In section III.B, we investigate the behavior of static two point correlations of the 
density fluctuations. Following the arguments presented in section II, non-equilibrium constraints are implemented by 
randomly initializing the nodes at the automaton edge s which are in contact with particle reservoirs. This procedure 
corresponds to assuming a local equilibrium form for the fluctuations at the boundaries. In the bulk, the fluctuations 
are seen to exhibit correlations that decay only algebraically over distances comparable to the size of the system. 
These long-range correlations are very sensitive to different implementations of the non-equilibrium constraints. This 
point is illustrated in section III.C, where the constraints arc imposed by direct manipulation of the microscopic 
configuration of the automaton. Simulations show that the qualitative features of the long-range correlations are 
modified, even though the macroscopic density profile is the same in both cases. 



A. Macroscopic dynamics: Diffusion equation. 

The starting point for the derivation of an equation for the time evolution of the average occupation of a given 
node is the equation that describes the microdynamics of the automaton 

ni(r + Ci,i+l) = n{N{r,t)) (26) 

where N{r,t) = '^^^^ ni{r , t) is the occupation number of node r at time t, and TZ{N) is a random fimction, 
which takes the value with probability (l — -j) or the value 1 with probability -j, where N is an integer 
between and 4. The Boltzmann equation is obtained by averaging Eq. (26) over the stochastic dynamics 

1 ^ 

/i(r,i+l) = ^ E /iCr-Ci.i), (27) 

3=0 

with the definition fi{r,t) = (nj(r, t)). The angular brackets denote the average over the stochastic dynamics. Sum- 
ming Eq. (27) over the index i, we obtain the equation governing the evolution of p(r, i) = {N{r,t)) = X]j=o 
the local density (or average node occupation number) 
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1 ^ 

+ = 4 E (28) 



i=0 



Note that this is an exact Boltzmann equation obtained without appeahng to the molecular chaos hypothesis. The 
stochastic dynamics contains intrinsically the factorization of the occupation densities through the random redistri- 
bution of the channel occupations [ p8[ . Equation (28) can be rewritten in the form 

p{r,t + l) - p{v,t) = \jz {e-^-^ -l)p{v,t). (29) 

1=0 

In the continuous time limit and expanding to second order in the gradients (hydrodynamic limit) Eq. (29) yields the 
diffusion equation 

|^p(r,i) - D\7lp{v,t), (30) 

with a diffusion constant D — ^. In the derivation of Eq. (30) we have made use of the symmetries of the lattice, 
namely 

3 3 

= 0; ^Cid = 2T, (31) 

i=0 i=0 

where 1 is the unit tensor. 

For the sake of comparison of theoretical results with simulation data, it is convenient to define reduced quantities, 
which depend only on one of the spatial directions, X, the direction along which the gradient is imposed, and on 
time. This is achieved by averaging over the remaining spatial direction, Y . Along this direction periodic boundary 
conditions are imposed by identifying nodes at y — Q and y ~ Ly. Hence, the reduced density per node is 



pi^'t) = T ^ P^""'^^- ^^^^ 

In Fourier space, one can easily write the time evolution 

5p{k,t) = Sp{k,0) exp{-Dk^t}, (33) 

27r 

k — -j— n, n = 0,1,2,... 

where 6p{k,t) is the Fourier transform of 6p{x,t) = p{x,t) — (a;) , the deviation of the actual density profile from 
the asymptotic stationary profile Ps{x). This stationary profile is a solution of the differential equation ^^Ps{x) = 0, 
and, depending on the boundary conditions imposed along X, it can be a homogeneous equilibrium profile Ps{x) — Peq, 
when p{0) = p{Lx) — Peq, or a non-equilibrium steady state, Psix) = Ps(0) + P x; (3 — BiLsl-J^^^ when 
p(0) ^ p(i,). 

The simulations of the automaton dynamics have been carried out under both equilibrium and non-equilibrium 
conditions. Figure 3(a) shows the time evolution of the reduced density profile in a (256 x 2048) automaton with 
periodic boundary conditions in both the X and Y directions. Initially the automaton exhibits a square density 
profile, in which all channels of all nodes are occupied for x < 26. This initial non-equilibrium profile evolves to a 
homogeneous equilibrium state characterized by a density per channel deq = 0.0976 [Note: The density per channel 
is just one fourth of the density per node: d{x, t) = ^^i^ ]. 

The evolution of the different Fourier components of the density profile p{k,t) is shown in Fig. 3(b). The amplitude 
of the mode fc = is constant, reflecting the global conservation of the number of particles. The remaining modes, 
with fc ^ 0, have an average exponential decay that agrees with expression (33). From this exponential decay we can 
measure the diffusion constant, which has an experimental value of D — 0.252 ± 0.002 p9|] , in excellent agreement 
with the value D = j, predicted by the theory. The deviation arises from the terms neglected in approximating a 
finite difference equation (Eq. (28)) with a differential equation in both space and time (Eq. (30)). 
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FIG. 3. Simulation data for a 256 x 2048 automaton with periodic boundary conditions in the X and Y directions, (a) 
Evolution of the reduced density profile along the X direction. Initially the automaton presents a square density profile with 
all nodes with a; < 26 completely occupied. Asymptotically the automaton evolves to an equilibrium state characterized by 
a constant density profile, (b) Time evolution of the different Fourier components of the reduced density profile. In this and 
subsequent figures time is given in automaton time step units, and space in units of lattice spacing. 

The results of simulation under non-equilibrium conditions in a somewhat smaller automaton (256 x 256) are shown 
in Fig. 4. The initial state of the automaton is a square profile where the nodes on the left half of the lattice are fully 
occupied. The evolution quickly renders the profile smooth and eventually leads to a linear density profile, which is 
the asymptotic stationary state for the implemented boundary conditions (d(x = 0) = 1, d(x — L^) — Q\ It is clear 
from Figs. 3(a) and 4 that the fluctuations in the 256 x 256 automaton are noticeably larger than in the 256 x 2048 
automaton, in accordance with the idea that the relative size of fluctuations decreases with increasing system size. 
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FIG. 4. Time evolution of the reduced density profile along the X direction in a 256 x 256 automaton under non-equilibrium 
conditions. Initially the system has a square profile with all the nodes for < a; < 128 fully occupied and all other nodes 
128 < a; < 255 empty. Asymptotically the system evolves to a steady state characterized by a linear density profile. Periodic 
boundary conditions are used for the Y direction. For X, the density at the boundaries is kept constant at the values 
d{x = 0) = 1 and d{x = 255) = 0. 



B. Mesoscopic dynamics: Fluctuations. 

We have derived an exact Boltzmann equation, Eq. (28), which describes the dynamics of the single particle 
distribution of the automaton. In order to characterize the fluctuations it is necessary to investigate the evolution of 
quantities associated with several particles. In particular, for two particles, we can formulate the equation 

n,;(r + Ci, i + 1) nj(r +Cj,t+l) 

= (l-<5(r,r')) E 'T,a;e{,}^we{.'}^{.'}5({s},Wr,i)})5({s'},{n(r',i)}) 

M.W} {s'},W} 

+ 6{r,r') J2 <y^<yJ ({s}^{a} S{{s},{n{r,t)}). (34) 

The summations are over possible configurations of a node, with the notation {s} = {s^; i = 0, 1, 2, 3}. The random 
function ^{^^j^^g.} is 1 or depending on whether the collision produces the configuration {a} from the initial 
configuration {s} or not. Upon averaging over the stochastic dynamics this quantity becomes A^gj_,^„j, a matrix 
which has a block-diagonal form, provided that the configurations {s},{(t} are grouped into equivalence classes 
according to the value of the conserved quantities. For the automaton considered here 

A{s}^{a} = — (5(^Sfc,^crfc), (35) 

'^{"J k k 

where u^^y is the number of elements in the equivalence class to which the configuration {s} belongs. Its value for 
the present case (see Appendix B) is 



Appendix B also contains a derivation of the result 



1 ^ 



4 

{a} k=0 
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which is needed, together with the Boolean nature of the variables ( af = ) to evaluate the average of Eq. (34) 
over the stochastic dynamics: 



ni(r + Ci,t + l) nj(r +Cj,t + 1) 
= (l-^(r,r)) l^(n,(r,i)nKr,t) 



k,l 



(5(r,r 



J% ^(nfe(r,i)) + ^ (l-<5y) J2{nk{r,t) (n,(r,i)-l)) 

k k,l 

1 (7V(r,i)7V(r ,t)) + 5(r,r')l (1-45,,) [(7V2(r, t)) - 4 (iV(r, i))] . 



(38) 



Summing over the indices i,j and making use of the equation for p{ic,t) = {N{r,t)) 
(Eq. (28)), we can derive the equation for the correlation function of the density fluctuations 

C{v,v';t) = ^N{v,t)Niv' ,t)) - (7V(r, i) ) (iV(r', i)). 



g-ci-V. g-Cj V- _ c{r,r;t) 



C(r,r ;i+l)-C(r,r',i) = ^ ^ 

ij 

+ i ^^"'^'■'^^^"'^'■'^^'(1 - 4%)[(7V^(r,<)) - 4p(r,i)]<5(r,r'). 



(39) 



This exact equation can be approximated in the hydrodynamic limit (small gradients) and in the continuous time 
limit by the following equation 



d_ 
dt 



C{v,v-t)^D (V2 + V^,) C{v,v';t) 



9n , 

— {{N\v,t)) 4p(r,t)) <5(r,r) 



(40) 



At this point, it is interesting to compare expression (40) with the theory of fluctuating hydrodynamics. Following the 
prescription given by Landau and Lifshitz |30| , we assume that dN{r,t) = N(r,t) — p{v,t) is a random variable 
obeying a stochastic equation, which is constructed by adding a noise term to the macroscopic diffusion equation that 
governs the evolution of p{r,t) = {N{r,t)) 



d_ 

di 



SN{r,t) = DV^dN{r,t) + V-g(r,t). 



(41) 



The term g(r, t) is a random particle flux assumed to be Gaussian white noise with the properties 

(g(r,t)) =0 

Jg(r,t) g(r',t')) ===T FFHir,t) 5(r,r') 5(1-1), (42) 

where Fp}i{v, t) is the amplitude of the noise. Given this specification for the noise term, we can derive the equation 
for the pair correlation function for the density fluctuations, C(r,r;<) = /dN{r,t)SN{r,t) 



d_ 

dt 



C(r,r;t) ^ D (V^ + VM C(r,r;0 + Vr • V^^ F^^^^ (r, t) 5(r, r 



(43) 



Knowing that, at equilibrium (i.e. for a constant density profile p(r, t) — p ), the correlation function for the density 
fluctuations is exactly (see Eq. (3.26) of Ref. [ ^ ) 



C"=«(r,r ) 



(4 - p) 5{r,r), 



(44) 



we conclude that the strength of the noise term must be 
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(45) 



in order that Eq. (43) give the correct answer for the equihbrium correlations. 

Away from equihbrium, we make use of the hypothesis that the noise term, which reflects the effect of dynamics at 
the microscopic scale, has a local equilibrium form 



FFH{v,t) = ^ p{r,t) (4 - p{r,t)). 



(46) 



With this hypothesis Eq. (43) becomes 
d 



dt 



C{v,v'-t) = D {Vl + V^,) C(r,r';i) 



jp{r,t) (4 - p{r,t)) <5(r,r') 



(47) 



We observe that this equation, which has been constructed in a phcnomcnological fashion, can be obtained from Eq. 
(40), which was derived directly from the microscopic dynamics, if we assume that the on-site correlations have a 
local equilibrium form 



3 



{N\r,t)) « {N'{r,t))^^ = - p^{r,t) + pir,t). 



(48) 



This is a good approximation in the thermodynamic limit, since the term neglected is of the order of one over size of 
the system. 

In order to study in detail the nature of these correlations, we focus on a {L^ + I) x Ly automaton in a non- 
equilibrium steady state characterized by an average stationary density profile, which is linear in the X direction 



Ps{r) = Ps{x = Q) + 13 X, 



(49) 



where /3 = ^ {p{x = L^) - p(x = 0)). 

In the pair correlation function for the density fluctuations we can distinguish two contributions: A local equilibrium 
term and a part containing the long-range correlations 



where 



C(r,r;i) = C^^^^(r,r) + C^^(r,r;t) 



C^^^^(r,r') = ip,(r) (4 - p,(r)) 6{r,r). 



(50) 



(51) 



Substituting expressions (50-51) in Eq. (39), we obtain 

C^>,r';i + 1)-C^>,r',i) = 1 E (e""' 



l) C^«(r,r';i) 



1 V- 



48 



(1 - 4Sij) C^«(r,r;t) S{v,r ) 
1 



p,(r) (4 - p,(r)) 



<5(r,r). 



(52) 



Using Eqs. (31) and (49) to simplify the last term, we have 

C^>,r ;t+l)-C^«(r,r',t) = 1 ^ (g-'^' V. e-'^J''^.' - l) C^>,r;t) 



+ ^ J^^-''''^' e-'^y^^' {1 - 4 5,,)[c^«(r,r;i)<5(r,r') 



48 



(53) 
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If we neglect the second term on the right-hand side of Eq. (53), which is equivalent to assuming that the on-site 
correlations have a local equihbrium form (i.e. {N^{r, t)) « (-/V^(r, i))^^ ), we find that this equation describes a 
random-walk in a four-dimensional cubic lattice (r, r ) with sources of constant strength along the surface r = r . 
In order to visualize the mechanism by which the long-range correlations are built up, consider the evolution of a 
system initially at local equilibrium ( C^^(r, r ]t = 0) = ), with a stationary density profile given by Eq. (49). 
On average, channel i of node r has an occupation /i(r) = i ps(r), for i — 0,1,2,3. After the propagation 
step, the average occupation per channel at node r is 



/o(r) 


1 

~ 4 


(Ps(r) 


- (3 Ax) 


/i(r) 


1 

~ 4 


Ps{r) 




/2(r) 


1 

~ 4 


(Ps(r) 


+ (3 Ax) 


/3(r) 


1 

~ 4 


Ps{r), 





where Ax is the lattice spacing in the X direction, which we have taken to be equal to one throughout the paper. 
The pair correlations are local and equal to {6ni{r) Snj{r)) = Sij /i(r) (1 — fi{r)), which yields 

C{r,r) = lp,(r) (4 - p,(r)) - i (54) 

After the collision step the particles are redistributed at random, which implies that, on average, the channel oc- 
cupation number is /i(r) = j /Os(r), for i = 0,1,2,3, as initially. Since the number of particles at each node 
is conserved in the collision step, the total correlations (Eq. (54)) do not change, except that the non-equilibrium 
contribution is now completely off-diagonal in the channels: 

C,,(r,r) = {dn,{r) Sn,{r)) = % /.(r) (1 - /,(r)) - (l-^,,)!/?^. (55) 

In C'ij (r,r) the diagonal terms (i — j) have a local equilibrium form and they remain local upon propagation. 

The off-diagonal terms { i ^ j ) give rise to the source term proportional to — ^ that appears in Eq. (53). Upon 
propagation, these off-diagonal correlations become non-local and perform a random walk (diffusion in the continuous 
limit) in the four dimensional lattice (r, r ). This accounts for the first term that appears on the right-hand side of 
Eq. (53). The only contribution unaccounted for is the second term in the r.h.s. of Eq. (53), which corresponds to 
long-range correlations created from pre-existing long-range correlations. It is small in the thermodynamic limit, and 
we shall neglect it henceforth. 

Finally, after neglecting the on-site contributions from the long-range term, and in the hydrodynamic and continuous 
time limits, we approximate Eq. (53) by 

|^C^«(r,r;0 ^ D {^l + V^,) C^«(r, r ; t) - ^-0'5{v,v). (56) 

In order to compare with the simulations, we define reduced quantities by averaging over The Y direction (orthogonal 
to the direction of the gradient) 

C{x,x-t) ^ ^(r^r';*)- (57) 

We impose periodic boundary conditions in the Y direction and Dirichlet boundary conditions in the X direction. 
As discussed in the previous section, the Dirichlet boundary conditions (i.e. that the correlations at the edges of 
the system have a local equilibrium form), correspond to the system being in contact with two reservoirs of different 
chemical potential. This corresponds to the paradigm of the automaton implemented in our simulations, where the 
configuration of all the nodes at a; = 0, are reinitialized randomly at every time step, with the constraint that 
the average densities at a; = 0, L^: be constant (although, in general, different from each other). The asymptotic 
solution of (56) with these boundary conditions is [ p2p^ ] 

C^^^^(x,x) = C^^^^(x,x) + C^^'^^{x,x), (58) 
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with 



^ -^Ly ps{x) (4 - ps{x)) 5{x,x (59) 

and 

C^'^^^{x,x) ^ -^hi „ x) 9{x ~ x) + x{L,~ x) e{x ~ x)] . (60) 

Simulations have been carried out in 11 x 11 automata with both equilibrium and non-equilibrium constraints. 
Figure 5 (a) shows the density fluctuation correlations in a closed system with periodic boundary conditions in both the 
X and Y directions. It can be seen that together with the local-equilibrium term, there is a long-range contribution, 
which is constant and negative, as predicted by Bussemaker et al ]l^ . This long-range contribution is well-known 
in classical statistical mechanics and appears solely as a consequence of global conservation of particles (see 
Appendix C). 

More physically relevant simulations are carried out in automata with periodic boundary conditions in the Y 
direction and in contact with particle reservoirs at the edges perpendicular to the X direction. Figure 5 (b) shows 
the density fluctuation correlations between the middle node and the other nodes C"^^(x,x = 5) in an equilibrium 
system with a homogeneous density profile. The average density per channel in this automaton d = 0.5. Note that 
the correlations are local and exactly equal (up to sampling errors) to the value Ly ^ 11 (see Eq. (59)). The average 
was performed over 10^ time steps. 




FIG. 5. Equilibrium correlations of the density fluctuations, C^''(x,2; — 5), in an 11 x 11 automaton, (a) Periodic 
boundary conditions in both X and Y directions; the average density is d = 0.5. (b) Periodic boundary conditions in the Y 
direction; for the X direction, the average density of the boundary nodes is kept fixed at the value d{x = 0) = d{x = 10) = 0.5. 

Figure 6 ^a) contains a plot of the correlations between fluctuations at nodes a; = 5 and nodes with arbitrary x, 
(jNESS^^^ X =5), for a system that is maintained in a non-equilibrium steady state by fixing the boundary densities 
at the values d{x = 0) = 1 and d{x = 10) = 0. Figure 6(b) is a plot of the long-range component C^^^^{x, x = 5) 
of the non-equilibrium correlations depicted in Fig. 6(a). The full line is simulation data, obtained by averaging over 
10^ time steps. The dotted line corresponds to the approximation given by fluctuating hydrodynamics, Eqs. (58-60). 
A similar comparison is shown in Fig. 6(c) for C^^^^ {x, x = 3). 

Note that the equation for C^^{r,r ;t) (Eq. (56)) contains information about the magnitude of the gradi- 
ent (/3), but not about its direction. This means that correlations are also long-ranged along the direction per- 
pendicular to the concentration gradient. This is illustrated in Fig. 6(d), where we have plotted C{y,y' = 
5) — C^^^^ {x,y,x ,5) for the same automaton. The correlations are different from those depicted in 
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Fig. 6(b) owing to the asymmetry introduced by the type of boundary conditions (fixed in the X direction, periodic 
in the Y direction). The theoretical curve is obtained by solving Eq. (56) by Fourier analysis. 




0.0 2.0 4.0 6.0 8.0 10.0 0.0 2.0 4.0 6.0 8.0 10.0 

Nocle(x) Node(y) 

FIG. 6. Simulation data and theoretical results for an 11x11 automaton in a non-equilibrium steady state with d{x = 0) = 1, 
d{x = 10) = 0. (a) Correlations of the density fluctuations, C^^^^{x,x = 5). (b) Long-range contribution C'"^^^ {x,x = 5). 
(c) Long range contribution C^^'^^(a;,x =3). (d) Long-range contribution C^'^^^ {y,y =5). 

C. Effect of boundary conditions. 

Unless we have access to its microscopic configuration, the only way to impose a non-equilibrium density profile in 
a diffusive system is to establish contact with two particle reservoirs of different chemical potential. Provided that 
the diffusion coefficient is independent of density, a system under these constraints will display, on average, a linear 
density profile. Furthermore, as discussed in section II, the density fiuctuations exhibit both local-equilibrium and 
long-range correlations in the bulk, but only local-equilibrium correlations at the boundaries. 
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One of the advantages of investigating the type of problem under eonsideration with an automaton is that we are 
able to manipulate the configuration in each node at will. For instance, the system can be forced by injecting and 
withdrawing particles in such a way that it asymptotically displays an arbitrary stationary density profile. Obviously, 
this manipulation involves an important modification in the nature of the fluctuations, which are extremely sensitive 
to how the microscopic configuration of the automaton is altered, and the loss of detailed balance in those nodes 
where the external perturbation is acting. 

In order to illustrate the effect of different boundary conditions on the long-range fluctuations, we have carried out 
simulations in a diffusive LGA, where for < x < a stationary reduced density profile of the form 

Ps{x) = pi9{a-x) + {pi - 0x) 9{x - a)0{h - x) + p^Oix-h), (61) 

is established by creating particles at x = a, and annihilating particles at x = b at appropriate constant rates (see 
Fig. 7). This procedure is very different from the previous implementation, in which at every step all the nodes with 
x = a. h are re-initialized at random, in such a way as to maintain the respective average densities constant. As far 
as fiuctuations are concerned, this implied that the correlations between fluctuations in nodes at the boundaries and 
those in bulk nodes are systematically destroyed at every time step. As a consequence, long-range correlations in such 
a system are confined in the region a < x < b while the rest of the system { < x < a, b < x < ) exhibits only 
local-equilibrium fluctuations. 

For the implementation with the creation/annihilation procedure, in the continuous-time and hydrodynamic limit, 
the Landau equation governing the automaton evolution is 

^oix,t) = D^p{x,f) + ^g{x,t) + F,,t{x,t), (62) 
of, ox/ ox 

where the diffusion constant is D = j, the quantity g{x, t) is the familiar random particle flux with local equilibrium 
form, and Fext{x,t) is an external field representing the efl[ect of injecting and withdrawing particles at a; = a,b. 
The average of this external field is stationary and equal to 

F,xt{x) = {Fext{x,t)) = - D^psix) + f3D [5{x~a) - 5{x-b)\. (63) 

The external field has a non-vanishing fluctuating part, which can be assumed to be Gaussian white noise provided 
that we are only interested in calculating pair correlations. Since we are manipulating the system only locally at 
X = a,b, the fluctuations of the external field will also be local 

Fext{x, t) - Fextix)) {PextiV-, t ) - -Fext (?/)) ) 

= D5{t-t) [A„ 5{x - a)S{y - a) + 5{x - b)5{y - b)] , (64) 

where Aq (A?,) measures the intensity of particle creation (annihilation). The system is in contact with reservoirs 
at X = and x = of appropriate chemical potential so that the average density is constant and equal to pi 
and p3, respectively. 

Again, the correlation function of the steady state density fluctuations can be separated into two terms 

C^^^^''\x,a;') = C^^^^(x,x) + C^«^^(x,x'), (65) 
with a local equilibrium contribution, 

C^^^^(a;,a;') = \ Ly p,{x) (4 - p,{x)) 5{x,x'), (66) 
and a long-range contribution, which can be expressed as a double Fourier series 

oo oo 

C^^^^{x,x) = V V C„„ sin^x sin^x', (67) 

n=l m=l 



with 



11? + m? 



2X2 



2 Ly . mr . rmr 

-p sm— 2; sm-^^; 

2 Ja 



2 ./„ Lx Lx 



. mr . mr , mr , mr 

+ Ca sm 'Y^'^'- ^111 'Y^'^ ^ IT TT 
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where 



c„ = A„ - ^ /? (4 - 2pi) 
C6 = A, + J (4 - 2p3). 

Simulations have been carried in a 31 x 31 automaton with the parameters a = 10, b= 20, di = 0.75, ds = 0.25. 
For these values, Ca = Cb ~ c, for symmetry reasons. Figure 7 shows the density profile averaged over 10^ time-steps, 
which is indeed of the form given by Eq. (61). 



1 .00 




o.oo 



o.o 



10.0 20.0 

Node(x) 



30.0 



FIG. 7. Averaged density profile in a 31 x 31 automaton, obtained by creating particles at a = 10, and annihilating 
them at 6 = 20 at appropriate rates. The average densities at the boundaries are Icept at di = d{x = 0) = 0.75 and 
da = d{x = 30) = 0.25. 
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Node(x) 
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30.0 0.0 




10.0 20.0 30.0 
Node(x) 



FIG. 8. (a) Steady state correlations, C{x,x = 15), in a 31 x 31 automaton in the non-equilibrium state depicted in Fig. 
7. (b) Long-range contribution, C^"-^^{x,x' = 15). 
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In Fig. 8(a), the density fluctuation correlations C^^^^ {x,x ) axe plotted for x = 15. The continuous line in 
Fig. 8(b) connects the simulation data for the long-range correlations at steady state between an arbitrary point and 
the middle node C^^^^ {x, x — 15), obtained after averaging over 5 x 10^ time steps. The dashed line represents 
the approximate solution, given by expression (67) (the Fourier series has been truncated after n, m = 1000 ), for 
a value c = 9.1, which has been chosen so as to obtain the best fit to the simulation data. The agreement between 
simulations and the approximation given by Fluctuating Hydrodynamics is very good. There are several features 
which call for comments: First, the correlations extend beyond the central part of the system a < x < b, where there 
is a concentration gradient: They have a non-zero value for all < a; < L^- Second, the form of the correlations 
is qualitatively different from those obtained in a system with a linear density profile for a < x < b, maintained 
by contact with particle reservoirs at the extremes. In particular they are positive, non-monotonic as a function of 
distance and furthermore they are non-linear. Nonetheless, the fluctuations in the "bulk" (i.e. around the central 
node) are qualitatively similar in both situations (compare Figs. 6(a) and 8(a), and Figs. 6(b) and 8(b)). 

IV. SUMMARY. 

We have constructed a lattice gas automaton with a set of fully stochastic collision rules, giving rise to diffusive 
behavior in the hydrodynamic and continuous time limit. This model system exhibits spontaneous fluctuations similar 
to those found in actual fluids The statistical properties of these fluctuations are the main objective of the present 
study. In the simulations the automaton is kept at equilibrium or away from equilibrium by placing it in contact with 
particle reservoirs of equal or different chemical potentials, respectively. From a computational point of view, these 
conditions are implemented by randomly initializing the nodes at the boundaries of the system at every time step. 
Under equilibrium constraints, the correlations of the density fluctuations are localized on each node, a result which is 
consistent with the observation that in fluids away from critical points or hydrodynamic instabilities, pair correlations 
decay on a microscopic lengthscale. The local nature o f equilibrium correlations is the result of a balance between 
noise sources, which is lost as soon as the system is removed from equilibrium. It is thus seen that in a non-equilibrium 
steady state the density fluctuations are correlated not only locally in an equilibrium-like manner, but also over large 
distances. This implies that the non-equilibrium distribution function is not factorized, in contrast to the equilibrium 
case. Both the local-equilibrium and the long-range contributions to correlations have been measured in computer 
simulations. The results of simulations have been compared to theoretical expressions derived from the stochastic 
microdynamics of the automaton. As a flrst approximation to these exact expressions, one can derive expressions 
obtained by means of the phenomenological theory of fluctuating hydrodynamics, which are in excellent agreement 
with the results of simulations even in fairly small automata. 

For the automata analyzed in this paper where the microscopic evolution preserves the particle number and where a 
preferred spatial direction is established by the imposition of a density gradient, global detailed balance is absent as a 
consequence the imposition of non-equilibrium constraints at the boundaries of the system. Long-range correlations of 
the type described here appear generically in dynamical systems with conserved quantities, having large scale spatial 
anisotropy and lacking detailed balance. 

ACKNOWLEDGMENTS 

It is a pleasure to acknowledge Irwin Oppenheim for fruitful discussions. This work was supported by the EC under 
contract ERBCHBG-CT93-0404. J. P. B. acknowledges support by the Fonds National de la Recherche Scientifique 
(FNRS, Belgium). P. G. is grateful to the Japan Society for the Promotion of Science for their support. 



18 



APPENDIX A: 



The object of this appendix is to obtain the spectral decomposition of the operators -^k{x), and ^k(x)^, 
where the thermal diflusivity is k{x) = ki 0{—x) + K2 d{x) 9{L — x) + ki 6{x — L), with the condition that the 
eigenfunctions are a superposition of plane waves in the limit \x\ 00. This corresponds to the setup depicted in 
Fig. 1 with semi-infinite reservoirs. 

Consider the eigenvalue equation 



92 



Pk{x). 



(Al) 



We require that k(.t) Pk{x) and its first derivative be continuous at x = 0, L, and that the Pk{x) form an 
orthonormal set (in the Dirac sense, that is we require that jTl^dx k{x) Pk{x) P^'{x) = 5{k-k) ). The 
eigenfunctions can be written as 

e{-x) 

e{x) 6{L - x) 

0{x-L). (A2) 





1 

Kl 


dk 


cos 


kx 

—;= + Ok ^ 




sin 


kx 

\/«r 


+ 


1 

K2 


Ck 


cos 


kx 

—p= + dk 1 




sin 


kx 


+ 


1 

Kl 


Cfe 


cos 


k{x-L) , ^ 
, — + Jfe 




sin 



First, we implement the continuity conditions at x = 0, and obtain the equalities 



Ck ^ ak, 
bk — dk- 



The continuity conditions at x = L, yield the result 

kL 



efe 

fk 



Finally, the orthonormalization condition is 

dx k{x) Pk{x) Pfe' (.t) 



cos ■ , 
h= sin -7^ 



K2 sm 



kL 



cos 



kL 



Ck 

dk 



I 

O — ( 



Kibl + el + Kifl) 5{k-k) 



+ {aka-k - Kibkb-k + ekS-k - Ki/fe/-fe) <5(fc + fc )| = 5{k-k) 
Combining Eqs. (A3-A5), we obtain 

bk = Oi\k\ dk, 

with 



1_ ]_ 

K\ 1 



^2kL 



cos 



2 fcL 



«2 sin2 fei 



(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



We are mainly interested in the limit e = , / — 



; then 



a 



\k\ 



1 1 

Kl K2 1 



^2 fcL 

\/k2 



,2 _fei^ ' 



(A8) 



when e ^ 0. Note that if ^ n it for any n d Z then all coefficients a^, fk are at most of order e. On 



the other hand, if 



kL 



ni: + T], with 77 <C 1 for some n G Z , then all coefficients are at most of order e 



except afc, which has a singular behavior. 
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2 V«2 



TT (2e2 + ?72) - ^ ?y (2e2 + r?)3 
For e — > Eq. (A9) is an asymmetric representation of a delta function. That is, 

4 ~ ^/K^ 6{r}), 

when e ^ 0. 

Thus, to lowest order in e, the eigenfunctions of the operator are 

Pk{x) = ^ cos ^(a;) ^(L - a;), 

K2 jKi 



with 



kL 



— mr) 



) = f E ^(^ - 



71" 



)• 



(A9) 



(AlO) 



(All) 



(A12) 



For the second operator, the eigenvalue equation is 



'dx'^^^^'dx ^^^^^ = - fc^ Qk{x), 



(A13) 



with the conditions that Qk{x) and i^{x)^ Qk{x) are continuous at x = 0,L. We further require that the set of 
eigenfunctions be orthonormal in the Dirac sense dx Qk{x)Q^.' (x) — 5{k — k). 

The solution of Eq. (A13) can be written as 



Qk{x) 



+ 
+ 



Qk cos 



Ck COS 



Cfe COS 



kx 




sin 


kx 








kx 


dk 


sin 


kx 








k{x - L) 


+ - 
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sm - 



0i-x) 

0{x) e{L - x) 
e{x-L). 



Using the continuity and orthonormalization conditions, we have 



(A14) 
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efe \ 
fk J 
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cos 


kL 

sin^ 




= aik\ 
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4i 
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1 + 


C0S2-^ 
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C0S2-^ 



T= sin^ 
COS ■ ''^ 



Ck 

dk 



^ sin2-^ 



+ «2 sin^-^ 

K2 Kl \/K2 



(A15) 



In the limit e 



/ «2 ^ we have 



with 



Qk{x) w — sin-^ ^(a;)6'(L-a;), 

K2 \/K2 



a2 _ «2 



E '^(^ - 



(A16) 



(A17) 



when e 



0. 
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APPENDIX B: 



Here we compute the values of i^^gj and of J2{a}'^i ^^^{«}^{<^} fo^' t^^^ automaton described in Section III. 
We first evaluate the number of elements in the equivalence class of configurations with J2k particles: 

{a} k k 



{a} ^'^ I fc=0 J 



. ^277 3 



<7fc=0 



— / rfa; TT exp{i (sfc - (Tfe) x} 
da; expji E'Sfes:! (l + e"")** 



E 



« / 27r 



y" da; expji (Esfe - O^^l = E(1) '^(^'E'*'=) ^ ( Efe 



Thus, the number of configurations with a given number of particles is 



4 

EfcSfc 



(Bl) 



i.e., the number different ways of placing particles in the four different channels. We now evaluate the 

numerator of the second expression, which is the number of configurations with channel i occupied 

E o-i '^(E^'='E^'=) " 2~ / expii E^fc 2:1 (l + e"'^)^ HE e-'""*^ 

{ct} fe fe I fe J ki^iai=(S 



= h ^^pf = (e.'.-i 

(B2) 

This result could have been obtained straightforwardly by arguing that the quantity calculated is equal to the number 
of different ways of placing the remaining Efe ■Sfe — 1 particles in the 3 available channels, given that channel i is 
already occupied by one particle. 
Combining (Bl) and (B2), we have 



E{<x}'^» ^(Efe-sfe.Efe^fe) _ Y^k 



^"^^^^^^^^^ " Ew ^i^kSk,Ek-k) 



Sk 



4 ' 



(B3) 



which just states the fact that after a collision, the average occupation of a channel is proportional to the number of 
particles participating in the collision. 



APPENDIX C: 

In a closed system at equilibrium long-range correlations appear rather trivially as a consequence of global conser- 
vation laws. For the automata under consideration, it has been verified (see Fig. 5(b)) that equilibrium correlations 
of the density fluctuations exhibit a long-range contribution, which is constant, and whose integrated value is equal 
in magnitude and opposite in sign to the local equilibrium term, so that the total number of particles be conserved. 

Consider an automaton whose lattice C contains V nodes. There arc b cqiiivalcnt channels per node. The particle 
density (number of particles per node) is p. We consider periodic boundary conditions, so that the system is closed 
(i. e. the total number of particles N = pV is conserved). Let {n{-)} = {nj(r); r e £, i = 0,1. ..(6— 1)} denote 
the automaton configuration. The equilibrium probability of having the configuration {«(•)} is 
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p[{nm = n n pkw) M EE - A ' (CI) 

r i \ r i / 

with p(nj(r) = l) = 6, p{ni{r)=0) = 1 — 6. The parameter ^ is determined by normalization: 
E PiM-)}] = E n n ^^('^^W) i: r expLa;(^^n,(r)) - ivl 

{n()} {«(•)} r i •'O I \ r i / J 

= i- /'^ e-^ n n E P(^iW)e^^-(^) 

° r * n»(r)=0 

1 /■27r 

= io "^"^ n n [^^" + (1 - ^)] 

= ^ dx e-"^ r^e^^ + (1 - 0)1 

Jo 



TV 



The last equality yields 9 as an implicit function of the automaton parameters V, b and of the occupation A''. 
In a similar way, the average occupation per channel is 

(n,(r)) = J2 n.ir)P[{n{-)}] = ^ T e— ^ [^g- + (1 _ ^)] ^""^ ^ e^x ^ P (^3) 

and the two particle distributions, for (r, i) ^ (r are 

;n,(r)n, (r')) = ^ n,(r) n,(r') P [{«(•)}] 

{"(•)} 

= ^ r dx e--^ [0e- + (1 - 6)] ^2 e^- = f (C4) 

ZTT V — 1 

By contrast, in an automaton which is maintained at equilibrium by contact with particle reservoirs, the total 
number of particles is conserved only on average, and the equilibrium distribution is 



p[{nm = n n "^(ee - A > 

V i \ r i / 



(C5) 



with p(ni{Y = 1)) = |, p(nj(r = 0)) — 1 — f , which guarantees normalization. The average occupation per 
channel is also {ni{v)) = ^, but correlations differ from the previous case (compare Figs. 5(a) and 5(b)). For 

(r,i) 7^ (r ,j), the two particle distribution factorizes (ni{r)nj {r )^ = (ni(r)) ^nj(r )^ = (|)^, which implies 
that the correlations are strictly short-ranged, in a system in contact with particle reservoirs. 
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